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Abstract 

This paper examines sparse grid quadrature on weighted tensor prod- 
ucts (wtp) of reproducing kernel Hilbert spaces on products of the 
unit sphere § 2 . We describe a WTP quadrature algorithm based on 
an algorithm of Hegland [TO] , and also formulate a version of Wasil- 
kowski and Wozniakowski's wtp algorithm [21j . here called the WW 
algorithm. We prove that our algorithm is optimal and therefore lower 
in cost than the WW algorithm, and therefore both algorithms have the 
optimal asymptotic rate of convergence given by Theorem 3 of Wasil- 
kowski and Wozniakowski [21]. Even so, the initial rate of convergence 
can be very slow, if the dimension weights decay slowly enough. 



1 Introduction 

This paper examines sparse grid quadrature on weighted tensor products of 
reproducing kernel Hilbert spaces (rkhs) on the unit sphere § 2 C M 3 . As per 
our previous paper on sparse grid quadrature on the torus [TT] , the empirical 
rates of convergence of the quadrature rules constructed here are compared 
to the theory of Wasilkowski and Wozniakowski [21] . 

The setting is the same as that used by Kuo and Sloan [15J, and Hesse, 
Kuo and Sloan [12] to examine quasi-Monte Carlo (qmc) quadrature on 
products of the sphere S 2 , but here we examine quadrature with arbitrary 
weights. 

As noted in our previous paper [TT], rates of convergence and criteria for 
strong tractability of quadrature with arbitrary weights are known in the case 
of weighted Korobov spaces on the unit torus [T3j [18] . As far as we know, 
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this paper is the first to examine the analogous questions for quadrature with 
arbitrary weights on the corresponding spaces on products of spheres. 

The algorithms we examine are an adaptation of the algorithm used in 
our previous paper [TUl E] , and an adaptation of the wtp algorithm of Wasil- 
kowski and Wozniakowski |21j. We examine these algorithms theoretically, 
giving bounds for the asymptotic convergence rate of quadrature error in the 
worst case. We also examine the algorithms empirically, via a small number 
of numerical examples. 

Arbitrary weight quadrature on products of spheres is interesting, not just 
for purely theoretical reasons, such as comparison with the results for equal 
weight quadrature, but also for practical reasons. Integration over products 
of the unit sphere is equivalent to multiple integration over the unit sphere. 
Such multiple integrals can be approximated in a number of ways, including 
Monte Carlo methods. 

Applications of tensor product spaces on spheres and approximate in- 
tegration over products of spheres include quantum mechanics [23], and 
transport and multiple scattering problems in various topic areas, including 
acoustics [TT], optical scattering problems [DEUE], and neutron transport 
problems [20] . 

One prototypical problem to be solved is scattering by a sequence of 
spheres. This can be modelled using a multiple integral of a function on 
the product of the spheres. The decay in the weights of successive spheres 
could model the decreasing influence of scattering on each successive sphere, 
as opposed to just cutting off the calculation after an arbitrary number of 
scatterings. 

The remainder of this paper is organized as follows. Section 2 describes 
our weighted tensor product space setting in detail. Section 3 describes the 
optimization problem involved in dimension adaptive sparse grid quadrature. 
Section 4 introduces the dimension adaptive (da) algorithm and shows that 
it is optimal in a certain sense. Section 5 analyses a version of the wtp 
algorithm of Wasilkowski and Wozniakowski, and compares its theoretical 
rate of convergence with that the da algorithm. Section 6 contains numerical 
results, comparing the two algorithms, and showing how the da algorithm 
performs as the dimension is increased. 

2 Setting 

The general setting used in our previous paper [TTJ applies equally well here. 
We repeat it here for the sake of completeness. 

Let T> C W +1 be a compact manifold with probability measure \x. It 
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follows that the constant function 1, with l(x) = 1 for all x G T>, is integrable 
and J v l(x) dfi(x) = 1. Then let if be a Hilbert space of functions / : T> — Y R, 
with inner product (•, •)#, and kernel K, with the following properties. 

1. Every / G H is integrable and 

/ = (1) 

Jo 

2. For every x E T>, the function fc^ G if, given by k x (y) := K(x,y), 
satisfies 

f(x) = {k x J) H , for all f E H. (2) 

We recognize if as a reproducing kernel Hilbert space (rkhs). In this frame- 
work, quadrature rules Q, defined by 

n 

Q(f) :=J2 w tf&) 
»=i 

are continuous linear functionals and Q(f) = (q, f)u with q = XT=i w ikxi- 

We will assume that the quadrature points Xi are given. An optimal 
choice of weights Wi minimizes the worst case quadrature error e(q), which 
is given by the norm ||1 — q\\ H . The optimal q* is thus defined as 

q* := argmin 9 {||l -q\\ H \ q G spaa{A: a . 1 , . . . , k Xn }} . 

The weights of an optimal quadrature rule are thus obtained by solving a 
linear system of equations with a matrix whose elements are the values of 
the reproducing kernel K(xi,Xj) = (k Xi , k Xj )n- The right-hand side of these 
equations is a vector with elements all equal to one. 

We now describe our more specific reproducing kernel Hilbert space % 
of functions on V. The space H satisfies fl2]), but as well as flTJ, it is also 
assumed to satisfy (the more specific) 

/ f(x)diM{x) = (l,f) H = 0. 
Jv 

We now extend "H into the space H 1 , which consists of all functions of the 
form g — al + f, where a G R, and / G H with the norm || ■ \\ n -y defined by 

hWu-y = l«| 2 + - WfWn- 
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It is easily verified that is an rkhs with reproducing kernel 

IC y (x,y) = 1 +iK(x,y), 

where K is the reproducing kernel of "H. 

For functions on the domain T) d we consider the tensor product space 
Hd '■= &>t=i H- lk where 1 ^ 71 ■ • ■ ^ 7d ^ 0. This is an rkhs of functions 
on V d with reproducing kernel K,d{x,y) := rifc=i(l + Ik fc{xk, Uk)) where 
Xk, Uk G T> are the components of x, y G V d . Moreover, since Tid is an RKHS, 
one has 

/ f(x) dn d (x) = (1, f) Hd , 

where /i^ is the product measure, (•, -)u d is the scalar product on the tensor 
product space Hd, and 1 is the constant function on T> d with value 1. 

The specific setting for this paper is that of Kuo and Sloan [15] . with 
s := 2. We now describe this setting. We take our domain T> to be the 
unit sphere S> 2 := {x G M 3 | x\ + x\ + x\ = 1}, and consider the real space 
L 2 (S 2 ). We use a basis of real spherical harmonics Y^ m (x), £ = 0, . . . , 00, 
m = —£, ...,£. For / G L 2 (S 2 ), we expand / in the Fourier series 

00 e 

f(x) = fo,o + ^2^2 ft,mYl, m (x). 
1=1 m =-t 

For a positive weight 7, we define the RKHS 

<J :={/:§ 2 ^M| H/IU) <oo}, 

1,7 

where 

00 1 

(/,ff) H W := /o,0 00,0 + 7 _1 (^ + 1 )) r /f,m3£,m- 

Kuo and Sloan [T5] show that the reproducing kernel of W[l is 

1 + jAr^x ■ y), where for z G [—1, 1], 

where Pi is the Legendre polynomial of degree I. Convergence of A r needs 
r > 3/2. 



Hi(x,y) := 
A r (z) := 
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For 7 := (7^1, . . . , 7d,d), we now define the tensor product space 

fc=i 

This is a weighted rkhs on (§ 2 ) d , with reproducing kernel 

fc=l 

Kuo and Sloan [15] studied equal weight (qmc) quadrature on the space 
H^, and found that it is strongly tractable if and only if Ylt=i Jd,k < 00 
as d — > 00. Hesse, Kuo and Sloan [T2] constructed sequences of QMC rules 
on this space, and proved that their worst case error converges at least as 
quickly as the Monte Carlo error rate of 0(n~ l l 2 ), where n is the cost of the 
quadrature rule in terms of the number of points. 

The work of Hickernell and Wozniakowski [13], and Sloan and Woznia- 
kowski [18J, on the weighted Korobov space of periodic functions on the 
unit cube, and the work of Wasilkowski and Wozniakowski [2T] on WTP 
quadrature on non-periodic functions on the unit cube, suggests bounds on 
the worst case error for our case of quadrature with arbitrary weights on the 
space H^. We might expect quadrature with arbitrary weights on this space 

to be strongly tractable if and only if Ylk=i ld,k < 00 as d — » 00, and, in the 
case of exponentially decreasing weights, as studied here, we might therefore 
expect the optimal worst-case error to have an upper bound of order 0(n~ r ), 
and a lower bound of order Q(n~ r ), for r > 3/2. 

Our analysis below shows that these expectations for quadrature with 
arbitrary weights on H^, can be met, and specifically, that our algorithm 
satisfies the upper bound suggested here. 

3 Optimization Problem 

We first describe the optimization problem in the general rkhs setting, as 
given in Section [2j 

We assume that a sequence of quadrature points x\,x 2 , ■ ■ ■ G P, and a 
sequence of positive integers n < n\ < . . . are given and are the same for 
all spaces "H 7 . The quadrature rules for "H 7 are then defined as some element 
of Vj' := span{/c 7 i; . . . , k^ n } C W 1 . We denote the optimal rule in by 
qj. Now define the pair- wise orthogonal spaces Uj by Uq := Vq , and by the 
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orthogonal decomposition V? +l = © Uj +l . Using the fact that the q- are 
optimal, one can see that 

P — P - P G TP 
and 5 := q$ <E Uq = Vq. Note that one has 

This is the fundamental reason why one needs the admissibility condition 
discussed in later this section. 

We use the notation J := N d , treating elements of J as indices, with 
a partial order such that for i,j G J, i ^ j if and only if ih ^ jh for all 
components. 

For a index i G J, let \.i denote the down-set of i, defined by \.i := {j E 
J I j ^ p. 13]. Subsets of J are partially ordered by set inclusion. For 
a subset / C JJ, let \,I denote the down-set of /, defined by \,I :— {J i€l li- 
Then J, J is the smallest set J D I such that if i G J and j ^ i then j G J. 
Thus HI = II. 

A sparse grid quadrature rule is then of the form 

je/ fc=i 

for some index set I. From the orthogonal decomposition V^' = @\ =l one 
derives the multidimensional orthogonal decomposition 

jell k=i 

One can then show that an optimal q G Vj is obtained as 

je4./ fe=i 

Thus both Vi and g| are obtained in terms of the down-set II, effectively 
restricting the choice of the set I to index sets which are also down-sets. 

This leads us to defining the concept of an optimal index set. An optimal 
index set is one which minimizes the error for a given cost, or minimizes the 
cost for a given error. Here, the cost is the number of quadrature points, 
which is the dimension of Vj. 
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We now make this definition precise. We first define vfj := dimUj*'* and 
:= 5] d ' k . For the remainder of this section, we use e G (0, 1) to denote the 

Jk Jk \ ' / 

required upper bound on quadrature error. The optimization problem then 
uses the following definitions. 



Definition 1. For index j G J, define 

d d 
IK A - <g>€- PiHIAill 2 , O: Pj /vj. 



k=l k=l 

For subset I C J, define 
Also, define P := 1 — e 2 . 



i/ere ; i/ie fci/i component of the index j. 



Due to the properties of Vj and Ay, i/ and p satisfy 

v jt v(I) G N + , < Pi < 1, < p(I) < 1, p(J) = 1. (3) 

We now consider the following optimization problem, posed as a mini- 
mization problem on the variable I C J. 

Optimization Problem 1. 

Minimize v(\,I), subject to p(I) ^ P, 

for some < P < 1, where v and p satisfy (j3J). 

In other words, given a required upper bound e on the quadrature error, 
the problem is to find the subset I C J with the smallest cost v{\,I) = 
^2 jeli u ii satisfying the constraint 1 — ^2j eI Pj ^ e 2 . 

Optimization Problem [1] can have multiple solutions, since for H, I, J C J 
if J = \.H = \.I and both p(H) ^ P and P then both / and J are 

solutions to Optimization Problem [TJ The following problem breaks this tie. 
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Optimization Problem 2. 



Maximize p(I) subject to I solving Optimization Problem^ 

The solution of Optimization Problem [2] satisfies an admissibility condi- 
tion. 

Lemma 1. If I is a solution of Optimization Problem^ then 

I = 11. (4) 

Proof. Let J = where I is a solution of Optimization Problem [21 and 
therefore of Optimization Problem [TJ Then I C J and thus J satisfies the 
constraints of Optimization Problem [TJ since pi > 0. Therefore J is also 
a solution of Optimization Problem [TJ since = If I $= <A ^ 

follows from pj > that > and so I cannot be optimal. Therefore 
1 = J. □ 

In view of the admissibility condition, we reformulate Optimization Prob- 
lem |2] as: 

Optimization Problem 3. 

Minimize z/(J), subject to 1 = ^1, and p(I) ^ P, 
for some < P < 1, where v and p satisfy ([3]). 

As pointed out by (e.g.) Griebel and Knapek [9], some sparse grid prob- 
lems can be formulated and solved as knapsack problems. The resulting 
solution is optimal in terms of estimated profit for a given "weight" . 

We call Optimization Problem [3] a down- set- constrained binary knapsack 
problem. Each item in the knapsack is a product difference rule. The profit 
for each item is just the squared norm of a product difference rule, and this 
can be calculated precisely. The "weight" of each item is the number of extra 
points the product difference rule contributes to the overall quadrature rule, 
assuming that the admissibility condition applies. The relationships between 
Optimization Problem [3] and other more well-known knapsack problems are 
described in more detail Section HI 

4 Algorithm 

The dimension adaptive (da) algorithm to choose the set / for T> = S 1 , 
described in [11] is quite general, and applies equally well to the case here, 
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Algorithm 1: The dimension adaptive (da) algorithm. 
Data: error e, incremental rules Aj and their costs Vj for j £ J 
Result: e approximation q and index set / 
J:={0}; g:=A ; 
while ||1 — q\\ > e do 

j := argmaXjjrj | i ^ I and / U {i} is a down-set}; 

/ := / U {j}; g := q + Aj ; 



where T> is S 2 . For the sake of completeness, we repeat the algorithm here, 

with a slight change in notation, as Algorithm [TJ 

This is a greedy algorithm for Optimization Problem [3j 

Under certain conditions on v and p, Algorithm [1] solves Optimization 

Problem [3J To see this, consider monotonicity with respect to the lattice 

partial ordering of JJ. 

Definition 2. The function p £ is monotonically decreasing if i < j 
implies that pi ^ pj. If i < j implies that pi > pj, then p £ M.\_ is strictly 
decreasing. The definitions of "monotonically increasing" and "strictly in- 
creasing" are similar. 

Using Definition &\ the following theorem holds. 

Theorem 2. If p £ is strictly decreasing and v £ is monotonically 
increasing, then Algorithm^ yields a quadrature rule q and index set I such 
that I solves the down- set- constrained knapsack Optimization Problem^ for 
P = p(I) = ||l-gf. 

The proof of Theorem |5] presented below proceeds in these stages. 

1. We introduce a related binary knapsack problem, and show that if / 
is a solution of the binary knapsack problem, and / is also a down-set, 
then / is a solution of Optimization Problem |3j 

2. We define the efficiency r i; describe a greedy algorithm for the binary 
knapsack problem (Algorithm [2] below), and show that if the efficiency 
is strictly decreasing, then each set I produced by the greedy algorithm 
is a solution of the binary knapsack problem, and / is also a down-set, 
and is therefore also a solution of Optimization Problem [3j 

3. We show that if the efficiency is strictly decreasing, then Algorithm |2] 
produces the same sequence of sets as Algorithm [TJ 
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A binary (0/1) knapsack problem [5] related to the Optimization Prob- 
lem |3] is: 



Optimization Problem 4. 



Minimize u(I) 



subject to p(I) ^ P, 



for some < P < 1, where v and p satisfy (j3J). 

Usually a binary knapsack problem is posed as a maximization problem, 
where the selection is from a finite set of items. Here we have a minimization 
problem and a countably infinite set. A finite minimization problem can 
always be posed as an equivalent maximization problem [TBI P- 15]. In the 
case of Problem H] this cannot be done, because the quantity to be maximized 
(the sum of the costs of the elements not in the knapsack) would be infinite. 
Instead, we must deal directly with the minimization form. 

We now formulate a converse of Lemma [TJ 

Lemma 3. If I is a solution of the Optimization Problem^ and I also sat- 
isfies the admissibility condition 1 = ^1, then I is a solution of Optimization 
Problem^ 

Proof. I satisfies the admissibility condition \,I = I and consequently p(I) ^ 
P. It follows that I satisfies the constraints of Optimization Problem [3] and 
thus minimizes p under these constraints, i.e., is a solution of Optimization 



This justifies our calling Optimization Problem [3] a down-set-constrained 
knapsack problem. 

If, in Optimization Problem [TJ we identify each set / C J with its indicator 
function X G {0, 1} , where Ij = 1 if and only if % G J, we obtain a more 
usual formulation of the binary knapsack problem: 

Optimization Problem 5. 



for some < P < 1, where v and p satisfy (J3J). 

Solving the binary knapsack problem is hard in general, but for certain 
values of the constraint P, a greedy algorithm yields the solution. These 
values are exactly the values for which the solution of the binary knapsack 
problem equals the solution of the continuous knapsack problem, which uses 



Problem [3J 



□ 



Minimize 




subject to 
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the same objective function v as Optimization Problem El and relaxes the 
constraints Xj G {0, 1} to Zj G [0, 1]. Dantzig [5 J gives a graphical proof of this 
for the classical binary knapsack problem - the finite maximization problem. 
Martello and Toth [16, Theorem 2.1, p. 16] give an explicit solution for the 
continous problem, and a more formal proof. 

The greedy algorithm for Optimization Problem H] is based on the effi- 
ciency Tj := Pj/uj. The algorithms generates the initial values of an enumer- 
ation J® of J, t G N + , satisfying 

The algorithm recursively generates Im from I(t-i), until for some T the 
condition 

p(I(t-i)) < p < PihT)) 

holds, where 

'w -Lb' w - 

The greedy algorithm is therefore as follows. 



Algorithm 2: The greedy algorithm for Optimization Problem HI 
Data: error e, incremental rules Aj and their costs Uj for j G J 
Result: e approximation g and index set / 
J:={0}; g:=A ; 
while ||1 — q\\ > e do 

j := argmaXjjrj | i ^ /}; 

J:=/U{j}; g: = g + Aj ; 



This algorithm has the following properties. 

Lemma 4. For any < P < 1, Algorithm^ terminates for some t = T. 

For each t ^ 1, £/ie set generated by Algorithm^ Im is the solution to 
Optimization Problem^ for P = p(I^). 

Proof. The algorithm terminates because j^' is an enumeration of J and 
therefore Ylt^iPjM = 1' since Y^tLi = 1, but P < 1. 

When p(I(t)) = P the constraints of Optimization Problem H] are satisfied. 
Furthermore, as the method used the largest r^, the objective function v 
is minimised for Optimization Problem HI A more detailed proof can be 
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constructed along the lines of the proof of Theorem 2.1 of Martello and Toth 

□ 

The construction of the enumeration used in Algorithm [2] requires sorting 
an infinite sequence and is thus not feasible in general, but, in the case where 
p is strictly decreasing and v is monotonically increasing, the enumeration 
can be done recursively in finite time. 

Here and in the following we say that % is a minimal element of a subset 
of J if there are no elements j < i in that subset. The minimum is thus with 
respect to the lattice defined by the partial order in J. 

Lemma 5. If p is strictly decreasing and v is monotonically increasing, at 
each step t > 1 of Algorithm^ the index produced by the algorithm is a 
minimal element of the set In_u '■— JT\ I(t)- Also = 0, and therefore 1^ 
is a down-set. 

Proof. If p is strictly decreasing and v is monotonically increasing, then r is 
monotonically decreasing. By construction, r^t-i) ^ r^t), so the enumeration 
must have < j®. It follows that j^ 1 ' = 0. 

For t > 1, since is an enumeration of J, no element occurs twice, 
and so £ ^(t-iy ^ny later element j( t+s ) in the enumeration cannot 
be smaller than so is a minimal element of Ift-i)- Since all elements 
smaller than j'w occur earlier in the enumeration, we must have |j C I(t-i) U 
{j^}- Therefore, if I(t-i) is a down-set, then so is L t y Since I(t-i) = {0}, by 
induction, Im is always a down-set. □ 

Corollary 6. For each t ^ 1, the set generated by Algorithmic lit) the 
solution to Optimization Problem^ for P = p(Im). 

Proof. This is an immediate consequence of Lemmas [3j H] and [5j □ 

The set M( t ) of minimal elements of /£ is finite. One can thus find 
= j with largest r» in this set. This is how Algorithm [T] finds the index 
j := argmaXjjrj \ i £ I and / U {i} is a down-set}, even in the case where 
the efficiency r is not strictly decreasing. 

In the case where r is strictly decreasing, we have the following result. 

Lemma 7. If the efficiency r is strictly decreasing, then Algorithm^ produces 
the same sequence of sets I as Algorithm^ 

Proof. From the proof of Lemma we have that if r is strictly decreasing, 
then for t > 1, j^' as per Algorithm [2] is always the element of M(t) which 
maximizes r. This is exactly j := argmaXj{rj | i I and /U{z} is a down-set}, 
as per Algorithm [TJ For both algorithms, im = {0}. □ 
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All the pieces are now in place for the main proof of this Section. 

Proof of Theorem^ From Lemma |5] we see that if the efficiency r is strictly 
decreasing, then each set / produced by the greedy algorithm (Algorithm [2]) 
is a solution of Optimization Problem [3j From Lemma [7] we see that if 
the efficiency r is strictly decreasing, then Algorithm [2] produces the same 
sequence of sets / as Algorithm [TJ 

If p is strictly decreasing and v is monotonically increasing, then since 
r j = Vil v i then r is strictly decreasing. Therefore each set / in the sequence 
produced by Algorithm [1] is a solution of Optimization Problem |3j □ 

It remains to show how to construct the set of minimal elements of 1^. 
To do so, we define S(i), the forward neighbourhood of i G J [HI P- 71] as 

S(i) := {j G J | i < j and (i ^ i < j => i = i)} , 

that is, S(i) is the set of minimal elements of {j G J | i < j}. 

Let e be the standard basis of M J . To construct Mm, start with Mm = 
S(i^) = S(0) = {ei, . . . , ed}- Then given M( t _i) and i®, one obtains M( t ) = 

(V)\{< (t) })uS(i (t) ). 
Note that 

(M(t-i) \ i^}) U = (M (t _ 1} U \ 

As the minimal elements of l£L are either elements of M(t-i) (but not i^) or 
elements of S(i^) we see that this set is equal to Mw. 



5 Error bounds 

We will now describe a second variant of wtp quadrature, 
identical to the sequence of quadrature rules g( DA ) described in Section [3] 
above, except that the order in which the incremental rules are added to 
this second variant is essentially the order used by Wasilkowski and Woznia- 
kowski [211 Section 5]. This variant uses criteria similar to those used by 
Wasilkowski and Wozniakowski [21j Theorem 3] , but adapted to our setting. 
These criteria are 

hj - 9j-i|L« ^ VjCD 3 , for all j ^ 1. (5) 

1,7 

(corresponding to Wasilkowski and Wozniakowski [2T| (39)]), and 

(j + l)D jp ^l, forallj^l, (6) 
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(corresponding to Wasilkowski and Wozniakowski J2TJ (36)]), 
for some D G (0, 1) and some positive C and p. 



As a consequence of (jSj), we have 



a, 7 



nil* 

fc=i 



(fc) 



^ b(d,j), where 



1-S 



k=l 



Let (^d,jfc)j fc = 1, . . . , d, be a sequence of positive numbers. In contrast 
to Wasilkowski and Wozniakowski [211 Section 5], we do not stipulate that 
= 1. Define 



Mi) n- 



l-So 



Jk 



(7) 



k=l 



We therefore have b(d, j) / j) — > as 1 1 1 1 x — > oo. We order the in- 
cremental rules in order of non-decreasing b(d,j)/£(d,j) for each index j, 
creating an order on the indices j^ ww ^(/i) . We adjust fc) so that this 
order agrees with the lattice partial ordering of the indices. We now define 



r(WW) . 



N 



._ |j( ww )(i) ? . . . , j( ww )(iV)}, and define the quadrature rule 



.(WW) . 



, f/ (ww) 



To obtain a quadrature error of at most e G (0, 1), we set 



N(e,d) := {j | b(d,j)/adj) > (e/C^r,)) 



where 77 G (0, 1) and 

Ci(d,i7):= 
Finally, we define 



(WW) 



(WW) 

j'6Jjv(e,d) 

We can now present our version of Wasilkowski and Wozniakowski's main 
theorem on the error and cost of wtp quadrature [211 Theorem 3]. 
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Theorem 8. Let rj € (0,1). Assume that a sequence of quadrature points 
Xi, X21 ■ ■ ■ € S 2 , and a sequence of positive integers n < n\ < . . . are given 
such that the corresponding optimal weight quadrature rules qj := gj G H^j 
satisfy © and (jSJ) for some D e (0, 1) and some positive C and p. Then 
the quadrature rule q^ W ^ defined by (jHJ) has worst-case quadrature error 
) e, and its cost (in number of quadrature points) is bounded by 

™tO^«.) (-) 

where 



C{d,e) 



e J 

C,i Ut=2 (1 + C^ k 2 /e d , k 9(k, e)) /(A;, e)' 



(1 -DP)(1-L)2) P /(2(1-^)) 
/ n 2^ \ 1/(2(1-*))) 

log (^i 2 /(a,(i - 1? 2 )) 1 /^)) nt 2 (/(<, 6))6-vd-,) 



-i + 



5y |_ x J + > we mean max(0,x). 

Wasilkowski and Wozniakowski's proof, with s := 2, a := 1, applies di- 
rectly to our Theorem [HI once the change in £<j,i is taken into account. 

Corollary 1 of Wasilkowski and Wozniakowski [211 P- 434] presents a 
simpler bound for the cost of their wtp algorithm, and their simplification 
also applies here. 

Corollary 9. For every positive 5 there exists a positive c(d, S) such that the 
cost of the quadrature rule o e ^ W ' ) defined by (jHD is bounded by 



costt^'xcwi) gy 



For exponentially decreasing dimension weights 7^, Theorem 4 of Wasil- 
kowski and Wozniakowski [21] shows that the g( ww ) rules are strongly poly- 
nomial. 

Our sequence of rules g( DA ) is more efficient than q( ww \ in the sense that 

g (DA) ig 

based on the optimal solution of the corresponding down-set-con- 
strained continuous knapsack problem, as explained in Section |H As a direct 
consequence of Theorem [8] and Corollary [91 we therefore have the following 
result. 
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Theorem 10. Let rj G (0, 1). Assume that a sequence of quadrature points 
Xi,X2, ... € S 2 , and a sequence of positive integers n < n\ < . . . are given 
such that the corresponding optimal weight quadrature rules qj := g] G IHl'f] 
satisfy © and (J6j) for some D G (0, 1) and some positive C and p. Let 
I(t),QM A ^ be the index set and corresponding quadrature rule generated by 
iteration t of Algorithm^ based on the rules qj, for sufficiently small error 
e = e . 

Then the quadrature rule q^ A ^ has worst-case quadrature error e(gS ) = 
6(t) := a/1 — p(I(t)), and its cost v{I(t)) is bounded by 



u(I (t) )^C(d,e it) )(—) 
\ e (t)J 



where C(d,e( t )) is defined as per Theorem^ As a consequence, for every 
positive 5 there exists a positive c(d, 5) such that the cost of the quadrature 
rule q^ A ^ is bounded by 

i/(J w )< C (d,«5) (— Y + . 
\ e (t)J 



6 Numerical results 

With the estimates given by our analysis in hand, we now compare these to 
our numerical results. 

Since our underlying domain T> is S 2 rather than S 1 , we need to change 
some of the details of the algorithm in comparison to the algorithm used 
for the torus [TT]. Specifically, we need a sequence of rules on a single 
sphere, which yields "good enough" worst case quadrature error with op- 
timal weights. Our choice of points for our numerical examples is a sequence 
of point sets, consisting of unions of spherical designs with increasing num- 
bers of points, and non- decreasing strengths. 

For the unit sphere § 2 , a spherical design [7] of strength t and cardinality 
m is a set of m points X = {xi, . . . ,x m } C § 2 such that the equal weight 
quadrature rule 

j m 

Qx(p) ■= 

m z — ' 

h=l 

is exact for all spherical polynomials p of total degree at most t. 

One difference between the constructions for S 1 and S 2 is that the nesting 
of spherical designs is not efficient. The union of two spherical designs of 
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strengths t\ and ti is in general, a spherical design whose strength is the 
minimum of t\ and t%. In the case of our numerical examples, the first 
design of strength is a single point. The next design of strength 1 consists 
of two antipodal points, so nesting is possible in this case. After this, the 
resulting unions of spherical designs, in general, have strength no greater 
than 1. 

For the numerical examples, a combination of (approximate) extremal 
(E) and low cardinality (L) spherical designs are used, according to Tabled! 
These approximate spherical designs were all provided by Womersley jU [22] . 



Index j 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


Type 


L 


L 


E 


L 


E 


L 


E 


L 


E 


L 


E 


L 


Strength t 





1 


1 


3 


3 


7 


7 


15 


15 


31 


31 


63 


Cardinality m 


1 


2 


4 


8 


16 


32 


64 


129 


256 


513 


1024 


2049 



Table 1: Strength and cardinality of approximate spherical designs used with 
Algorithm [TJ in the numerical examples. 

If we let nij := \Xj\, the cardinality of Xj, and let tj be the strength 
of Xj, then, for the sequence of spherical designs chosen for our numerical 
examples, the extremal spherical designs have rrij = (tj + l) 2 , and the low 
cardinality spherical designs have rrij = (tj + l) 2 /2 or rrij = (tj + l) 2 /2 + 1, 
and in all cases tj ^ y/TrT] — 1. It is not yet known a sequence of spherical 
designs satisfying this lower bound on strength can be extended indefinitely, 
but there is rigorous computational proof for tj up to 100 [3j- Also, it is now 
known that an infinite sequence of spherical designs exists with the required 
asymptotic order of strength, that is, there is a sequence of spherical designs 
of cardinality rrij = 0(t^), [2] but the proof is not constructive and the 
corresponding implied constant is still unknown. 

We now turn to estimates for rules on a single sphere, in order to use them 
with Theorem [HJ On a single unit sphere, our sparse grid quadrature rule is 
an optimal weight rule qj = qJ.(Sj), based on an increasing union of spherical 
designs, Sj := IJi=o^- ^ s wors t case error is therefore smaller than that 
of the optimal weight rule qJ(Xj) based on Xj, the largest spherical design 
contained in the union, which is, in turn no greater than the worst case error 
of the equal weight rule q7'^ CiMC \Xj) based on Xj, with weights 1/|X,-|, 

e 2 fe)<e 2 (g7(A,))^ e 2 (g- (QMC) (^))- 
According to Hesse, Kuo and Sloan [121 Theorem 4], for 7 = 1, we have 
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the bound 

e 2 (g^ MC )(X,)) ^ctf r , 

and therefore 

e 2 (ql{X 3 )) <: c(^m]-l)- 2r ^ C im j r , 

for some C\ > 0. 

For general 7, as per Kuo and Sloan [T5], we have 

3 h=l i=l 
trij mj 

= 7 ^2 E E ' x i>*) ^ (9) 

3 h=i i=i 

For the sequence of spherical designs chosen for our numerical examples, 
we also have 2 J ^ rrij ^ 2 J ' + 1. We therefore have 

e 2 ( gj )^7C 2 2^, (10) 

for some C2 > 0. 
Recall that 

e 2 fe) = 1 - hjf m (r) 

1,7 

= 1 - II^-iIIhM - ll<?j - ^'-iIIhW 
= e 2 (<?;-i) - Iki -^-iIIhM > 

1,7 

for j ^ 1, since — qj-i is orthogonal to g^-i. Therefore 

e 2 fe-i) = e 2 ^) + ||^ - 9i-i|| 2 M . 

1.7 

Since e 2 (oj) ^ 0, using (fTOj) we obtain ||g,- — g,_i|| 2 (r) ^ 7C2 2~ r: '~ 1 This, in 

' 1,7 
turn implies that 

1.7 

where C 3 = 2" r C 2 . 

All the approximate spherical designs listed in Table [I] have one point in 
common, the "north pole" (0, 0, 1). Therefore, in our numerical examples, the 
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number of points rij = \Sj\ satisfies rij — 1 + Xli=o( mi ~ -0 = ~3 + ^2i=o mi - 
Since 2 l ^ < 2* + 1, we have 2 J+1 — j — 1 ^ rij ^ 2 J " +1 , and so m,- ^ 
rij ^ Trij+i, for j ^ 0. Also, our numerical examples obtain a value for C3 
of approximately 1.453. In view of (Q, and the preceding argument, criteria 
© and © hold with D = 2~ r / 2 , C = C 3 ~ 1.453 as above, and p = 2/r. 

Our numerical examples use r = 3 and 7^ = g fc , for g = 0.1, 0.5, and 0.9, 
to see how our rules g( DA ) and g( ww ) behave as the decay of the dimension 
weights is varied. For the g( ww ) rules, we use £d,k '■= CD, with C and D 
defined as above. 

For the DA and WW weighted tensor product algorithms, each program 
run uses r = 3 ; g = 0.1 , 0.5, or 0.9; a particular dimension d, from d — 1 
to 16; a particular maximum 1-norm for indices, typically 20; and a partic- 
ular maximum number of points, up to 100 000. The numerical results are 
potentially affected by three problems. First, if 7 is close to zero, and the 
number of points is large, then the matrix used to compute the weights be- 
comes ill-conditioned, and the weights may become inaccurate. In this case, 
a least squares solution is used to obtain a best approximation to the weights. 
Second, if the current squared error is close to zero, and the squared norm for 
the current index is close to machine epsilon, then severe cancellation may 
occur. Third, the sequence of spherical designs used in our numerical exam- 
ples is finite, so it is quite possible that our algorithm generates an index 
corresponding to a spherical design which is not included in our finite set. In 
these last two cases, the calculation of the quadrature rule is terminated. 

Figure [1] displays the typical convergence behaviour of the DA and WW 
rules for the cases examined. The particular case shown is that of (S 2 ) 4 , 
r = 3 , 74^ = 0.5 fc . The number of points used varies from n = 1 to 100 000. 
The cost axis is horizontal and the error axis is vertical, to match the figures 
shown in the torus paper The curve in Figured] labelled "ww bound" is 
actually the minimum of the bounds given by Theorem [HI as the parameter 
7] is varied over a finite number of values between and 1. 

In general, the DA algorithm has a cost no greater than that of the WW 
algorithm. Both are bounded by the WW bound of Theorem [HI The WW cost 
bound itself has an asymptotic rate of convergence of (D(e~ p ) = 0(e~ 2 / r ) = 
0(e~ 2 / 3 ) for all of our cases. In other words, the asymptotic bound has 
quadrature error of order 0(n~ 3 / 2 ). Judging from Figure [lj the rates of 
convergence of both algorithms appear consistent with that of the bound, 
but the asymptotic rate is not achieved by either algorithm or by the bound 
itself, for the number of quadrature points displayed in the plot. 

For 7^fc = 0.1 fc , Figure [2] shows how the convergence rate of the error 
of the DA quadrature rules varies with dimension d, for d — 1 , 2, 4, 8, 
and 16. The curve for d — 1 appears consistent with the asymptotic error 
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J^Q— 5 I I I I I I I I I I I 

10° 10 1 10 2 10 3 10 4 10 5 
Cost (number of quadrature points, n) 

Figure 1: Error of DA and WW rules vs WW bound for (§ 2 ) 4 , r = 3 , 74^ = 
0.5 fc . 

rate e = 0{N~^/ 2 ). The cases d = 8 and d = 16 are almost indistinguishable 
on this Figure. This is an example of the convergence in dimension. 




Cost (number of quadrature points, n) 
Figure 2: Error of DA rules for (§ 2 ) d , d = 1, 2, 4, 8, 16 ; r = 3 , 7 djfc = 0.1 k . 
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Figure [3] shows the equivalent results for the DA quadrature rules for 7^ = 
0.9 fc . The curve for d — 1 again appears consistent with the asymptotic error 
rate e = 0(N~ 3 ^ 2 ), but as d increases to 16, the initial rate of convergence 
to zero of the error becomes much slower than that for 7^ = 0.1 fc . This 
behaviour is expected, given the WW bound. 




Cost (number of quadrature points, n) 
Figure 3: Error of DA rules for (§ 2 ) d , d = 1, 2, 4, 8, 16 ; r = 3 , ^ d)k = 0.9 k . 
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